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^ - Abstract 

r> . We calculate form- factor ratios between the semileptonic decays — > D~^i~i> and — > Dfi'u 

i with lattice QCD. These ratios are a key theoretical input in a new strategy to determine 

the fragmentation fractions of the neutral B decays, which are needed for measurements of 
BR(i3^ —7- fi^fi~). We use the MILC ensembles of gauge configurations with 2-1-1 fiavors of 
sea quarks at two lattice spacings of approximately 0.12 fm and 0.09 fm. We use the model- 
independent z parametrization to extrapolate our simulation results at small recoil toward maxi- 
mum recoil. Our results for the form-factor ratios are flf\M^)/flf'\M^) = 1.046(44)stat. (15)syst. 
and flf\M^)/flf\M^) = 1.054(47)stat. (17)syst.- In contrast to a QCD sum-rule calculation, no 
significant departure from [/-spin (d •(->■ s) symmetry is observed. 

PACS numbers: 12.38.Gc, 13.20.He 
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I. INTRODUCTION 



Recently there has been increasing interest in the rare decay — jj.^jj," which, as 
a flavor-changing neutral-current process, is forbidden at tree level in the standard model 
(SM). At the loop level, it can be mediated by weak bosons through penguin or box diagrams. 
With a nonperturbative (lattice-QCD) calculation of the bag parameter Bb^, the branching 
fraction has been predicted to be [1, 2] 

BR(E° ^ = 3.2(2) X 10"^ (1.1) 

Several new physics models would enhance the decay rate [3-7], and, hence, observation 
of this process could potentially reveal physics beyond the SM. Recently, several experi- 
ments [8-13] have published upper limits on this branching fraction, which we have compiled 
in Fig. 1. Moreover, CDF [11] reports an excess such that BR(SO fi+fi~) = 18^^] x 10"^ 
or a two-sided 90% confidence interval, 4.6 x 10~^ < BR(i?° — fi^jj,') < 39 x 10~^, lying 
above the SM prediction, Eq. (1.1). CMS and LHC6, however, set upper limits that restrict 
the CDF region. As statistics accumulate, especially at LHC6, a definitive measurement at 
the SM rate or higher seems likely soon. 

At a hadron collider, the extraction of BR(i?^ — )■ relies on normalization channels 

such as B^ — )■ J/tpK"^, B'^ — )■ K~^7t~ and B^ — )■ J/ip4' [14], through relations of the form 

BR(i?° ^ /iV) = ^ X)^ — (1.2) 

Is ^tifi 

where e and N are, respectively, the detector efficiencies and the numbers of events. The 
fragmentation fractions fq [q = u,d,s or A) denote the probability that a b quark hadronizes 
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FIG. 1. Comparison of the (most recent) measurements from CDF [8, 11], D0 [9], CMS [12], 
and LHC6 [10, 13] with the SM prediction [1, 2] shown as a vertical band. The filled bars show 
the measured bounds of the branching ratio with a 95% confidence. In the fourth bar, the inner 
box shows the two-sided 90% bound from CDF [11]. Two results from the LHC6 in 2011 are 
distinguished as "2011a" [10] and "2011b" [13]. 
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into a Bq meson or a Af, baryon. The fragmentation fractions fq may depend on the en- 
vironment, so they are best measured in situ in each experiment. Thus, improving the 
determination of the fragmentation ratio fs/fd will tighten the limits and increase the sig- 
nificance of measurements. 

The quantity fs/fd has generally been determined from semileptonic decays [15], an 
approach that LHCfe has newly refined [16]. Recently, Fleischer, Serra, and Tuning proposed 
two approaches based on measuring the ratio relative to nonleptonic decays BR(i?° — )■ 
L'+7r-)/BR(5° ^ D+K~) [17] or BR(50 -> Z)+7r-)/BR(50 D+vr") [18]. An important 
ingredient in both approaches is the approximate factorization of the nonleptonic decay 
amplitudes, which relies on the corrections to naive factorization of the light meson in the 
final state being small and calculable [19]. The D^K~ method is favored in this regard, 
because it receives contributions only from color-allowed tree-diagram-like topologies which 
yield smaller nonfactorizable effects [18]. 

The ratio BR(E° ^ L)+7r-)/BR(5° D+K') is related to fj fa by analogy with 
Eq. (1.2). Via factorization, the amplitudes for these nonleptonic processes can be expressed 
as a product of the light-meson decay constant and a semileptonic form factor for — )■ 
D(^s)^p. This leads to a way to measure fs/fd [17, 20]: 



^ = 0.0743 X — X 

fd Tbo 



^DK 

Ndk 



[1.3] 



where r denotes lifetimes, and the number 0.0743 is a product of ratios of well-known 
quantities such as the light-meson decay constants, Cabibbo-Kobayashi-Maskawa (CKM) 
matrix elements and kinematic factors. The factorization is parametrized by [17] 



Afp 



ai ' 



Ad) 



f!>'\Mj,) 



:i.4) 

:i.5) 



IS 



where a^'^^ is a factor accounting for the deviation from the naive factorization and fo^q"^ 
a form factor for the corresponding semileptonic decay. 

The hadronic method relies on theoretical inputs for A/^ and Afp. In the limit of exact 
[/-spin symmetry (namely the exchange of s and d quarks throughout the process), both 
reduce to 1. Fleischer, Serra, and Tuning expect the f/-spin breaking \Afa — 1| "to be at 
most a few percent" [18, 19]. Based on an estimate from QCD sum rules [21], they quote 
either Afp = 1-3 ± 0.1 [17] or Afp = 1-24 ± 0.08 [18], the latter of which LHC6 uses [20]. In 
either case, the biggest limitation is from the form-factor ratio Afp. 

A relation between fs/ fd and BR(i?° — )■ Df7i~)/BK{B^ — )■ D^n') is derived along similar 
lines [18]. In that case, the form-factor ratio becomes [/q*'*(M^)//q'^'*(M^)]^, i.e., with both 
numerator and denominator evaluated at = M^. 

In this paper, we calculate these two form-factor ratios using lattice QCD with 2+1 flavors 
of sea quarks. We use the same set of MILC ensembles of gauge configurations [22] and the 
same sequence of bootstrap copies for both of the and 5° processes, which reduces the 
statistical error by correctly accounting for correlations. We include the contributions of 
the first radially excited states in the fits of correlation functions to avoid the respective 
systematic errors. Such a treatment turns out to be necessary for calculations at nonzero 
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recoil. By fitting the correlation functions in a simultaneous and mutually constrained 
manner, we are able to extract the form factors at small recoil. We then extrapolate our 
lattice results to the continuum limit and to physical quark masses with the guide of chiral 
perturbation theory. Finally, we use the model-independent z parametrization [23] to extend 
the form factors toward large recoil. 
We finally arrive at the result 

^l,,^ = 1.046(44)(15), (1.6) 

where the first error is statistical and the second reflects the systematic errors added in 
quadrature. (Due to refinements in the analysis, Eq. (1.6) differs slightly from our prelimi- 
nary result [24].) We do not observe a large f/-spin breaking effect. Such a small difference 
between the and form factors is in accord, however, with recent lattice-QCD calcula- 
tions on lighter mesons like -D(s) — ti{K)Iv [25]. It is also in agreement with a result from 
heavy- meson chiral perturbation theory [26]. 

The factorization analysis of BR(i?° D^-k^) /'&K{B^ — D^tt^) is somewhat more 
complicated because of additional topologies in the decay B^ — D^ir". A similar form- 
factor ratio is needed and, simply by adjusting in the denominator, we find 



1.054(47)(17). (1.7) 



We discuss the implications of our results (1.6) and (1.7) in Sec. VII. Here we only note that 
both yield fragmentation-fraction ratios fs/fd in agreement with LHCfe's recent measurement 
via semileptonic methods [16]. 

This paper is organized as follows. In Sec. II, we summarize the formalism and our 
strategy to extract the form factors at nonzero recoil. We provide the simulation details 
in Sec. III. We describe the methodology that we apply to extract the form factors from 
the two- and three-point correlation functions with the given gauge configurations. This 
fitting procedure is crucial to our analysis. In Sec. IV we describe the chiral-continuum 
extrapolation using the corresponding chiral perturbation theory. In Sec. V, these results 
are then extrapolated to the region of small momentum transfer using a model-independent 
parametrization. We also compare here a related form factor, which we obtain as a by- 
product, with the experimental results. In Sec. VI, we account for the systematic errors 
that arise in our analysis and present a full error budget. Finally, in Sec. VII, we present our 
results, compare with previous results and discuss prospects and connections to current and 
future experiments. The Appendix specifies the functional form of the chiral extrapolation 
in detail. 



II. SEMILEPTONIC Di^^^lv FORM FACTORS FROM LATTICE QCD 

The hadronic matrix elements of the semileptonic decays Bi^g) D(^s)^^ can be parametrized 

by 



{D{p')\V^\B{p)) = U{q 



2^ 



A/f2 _ ]\/r2 



+ fo{q') ^ , ^ g", (2.1) 
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where q = p—p' is the momentum transfer and = cy'^b is the (continuum) vector current. 
Another parametrization uses velocity 4- vectors v = p/M instead of momentum p [27], 

<^'p;'i^''i^'^» = + .')" + ft-c-X" - "')". (2.2) 

^/MbMd +v a ; V A y , \ j 

where w = v ■ v' = (M^ + M|, — q^) /2MbMj:) describes the recoil of the process. The h± 
parametrization is convenient for lattice QCD, both for numerical simulation [28] and for 
matching lattice gauge theory to continuum QCD [29, 30]. 

The lattice-QCD calculation of h+ in the zero-recoil limit has been investigated using 
double ratios [28], 

_ {D\-c^%\B){B\h'c\D) _ 

^+ " {D\-c^^c\D){B\h%\B) ~ ^ 

with all states at rest. To proceed analogously at nonzero momentum, we introduce the 
following single ratios: 

^ {Dip)\cjb\B{0)) ^ h4w)-h.{w) 

{D{0)\n%\B{0)) 2h4l) ' ^ • ^ 

^ {Dip)\nb\B{0)) ^ h+{w)-h_{w) 

{D{p)\c-i%\B{0)) {w + l)h+{w) - {w - l)h_{w) ' 
^ {D{p)\-nc\D{0)) ^ 

~ {D{p)\c^^c\D{0)) l + w' ^ ' 

where the last follows from vector current conservation, h^~^^{w) = 0. 

We can write down the equations that manifest the relations between the ratios and the 
form factors 

1 + d- d , , 

hjw) ai 

ad, (2.8) 



h(w) ai Qi , , 

h+{l) h d- ^ • ^ 

In Eq. (2.4), we have a ratio a between matrix elements involving a final D meson with 
nonzero and zero spatial momentum. The purpose is to make use of the correlations in the 
uncertainties between the two. The form factor at zero recoil can be extracted precisely via 
i?+ [28], and we find that these ratios aid calculations at nonzero recoil in a similar way. 
With h±{w) in hand, one can obtain the form factors f+iq"^) and fo{q^), 

Uiq') = ^ [(1 + r)h+M - (1 - r)h.{w)] , (2.10) 



2\ 



1+r ^ ' 1 - r ^ ' 



(2.11) 



where r = Md/Mb and q^ = M% + MI-2wMbMd. Equations (2.8) and (2.9) both contain 
the factor /i+(l), so we write 

fo{q') = K{l)h{w{q')). (2.12) 
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TABLE I. Parameters of the MILC asqtad ensembles of configurations and the valence quarks used 
in this analysis. 

Ensemble a (fm) ami/arrih -/Vconfs csw amx{B — > D) amx{Bs — )■ Dg) 

C020 0.12 0.020/0.050 2052 0.1259 0.0918 1.525 0.020 0.0349 

C007 w 0.12 0.007/0.050 2110 0.1254 0.0901 1.530 0.007 0.0349 

F0124 w 0.09 0.00124/0.031 1996 0.1277 0.0982 1.473 0.0124 0.0261 

F0062 « 0.09 0.0062/0.031 1931 0.1276 0.0979 1.476 0.0062 0.0261 



In the formulas until now, we have not specified the spectator mass, so they apply to both 
the B ^ D and Bg — t- processes. With this notation, the desired ratio of the form factors 
is then 

where the first factor is obtained from the ratios R^~^^ and R^"'^^'' and the last from the 
expressions in Eqs. (2.8) and (2.9). 

On the lattice, we define a vector current = Zy^Zy^'^^l^'^b [28, 30], where the factors 
Z-iri normalize the flavor charge. The matching between the lattice and the continuum 
physics can be bridged by the relation = pvt^V^, where = Zyi Zy4 / Zyi Zy^^- The 
normalization factors Zy4^ cancel in the ratios in Eqs. (2.3)-(2.6). The factor py4, has been 
verified to be very close to 1 with one- loop perturbation theory with unimproved gluons [30] . 
Calculations of pyi^ with improved gluons (as used here; cf. Sec. Ill) are in progress. Given 
the ratio structure in Eq. (2.13), it is clear that the (small) contributions from pyt^ — 1 should 
largely cancel. Thus, in this analysis, we take pyt^ = 1 and estimate the uncertainty from 
this choice in Sec. VI. 



III. SIMULATIONS AND FITTING METHODOLOGY 
A. Data Setup and the Lattice Simulations 

Our calculation uses four ensembles of MILC's (2-|-l)-fiavor asqtad configurations [22] at 
two lattice spacings, a ~ 0.12 fm, 0.09 fm, which we refer to as the "coarse" and "fine" 
lattices, respectively. The configurations were generated with an O(a^) Symanzik improved 
gauge action [31-34]. The coarse (fine) ensembles used here have a lattice size of 20'^ x 64 
(28'^ X 96), so in both cases the spatial size is L 2.4 fm. The four ensembles have 
different sea-quark masses, so, for the sake of convenience, we label them C020, C007, F0062, 
and F0124. Details on the parameters that we use in the simulations are summarized in 
Table I. The strange and light sea quarks are simulated using the asqtad-improved staggered 
action [35-39]. The asqtad action is also used for our strange and light valence quarks. The 
heavy charm and bottom quarks are simulated using the Sheikholeslami-Wohlert (SW) clover 
action [40] with the Fermilab interpretation [41]. We simulate the B ^ D and Bg — >■ Dg 
decays on the same ensembles, so that correlations reduce the statistical uncertainty in the 
ratios. For the B D decay, the valence light-quark mass is taken to be the same as the 
sea-quark mass, i.e., we stick to "full QCD" data with nix = ''^z, while for the Bg — Dg 
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process, we set the valence strange-quark mass to be close to its physical value, rrix = mg. 
The charm and bottom quarks in our calculation are tuned to their physical values up to 
a tuning uncertainty. Columns 2-9 in Table I list, respectively, the approximate lattice 
spacings, light /strange sea quark masses, number of configurations, the hopping parameter 
/tb(c), the coefficient for the clover term csw? and the light valence-quark masses used in the 
B ^ D and — j- Dg simulations. The quark masses here are all in lattice units. 

We obtain the matrix elements appearing in Sec. II from the following three-point corre- 
lation functions: 

C3^^'*^(0,t,T;p) = (OD(0,0)^,z7'^vl/,(t,^)ot^(T,a.)) e^^^ (3.1) 
x,y 

C3^^^^(0,t,T;p) = {On{Q,0)^cil'^c{t,y)OUT,x)) e^^\ (3.2) 

x,y 

C3^^'^(0,t,T;0) = Y ((^B(0,0)W^fe(t,2/)(^UT,x)), (3.3) 
x,y 

where the sum over x sets the B meson at rest, and the sum over y selects the final-state D- 
meson momentum. The final D meson is simulated with several small spatial momenta which 
are the lowest possible values for the finite spatial volumes: p = 27r(0, 0, 0)/L, 27r(l, 0, 0)/L, 
27r(l, 1, 0)/L, 27r(l, 1, 1)/L, 27r(2, 0, 0)/L, and permutations. To increase statistics, data 
are generated at four different source times, spaced evenly along the temporal extent of 
the lattice. The zero-momentum correlation functions for D D and B ^ B serve as 
normalization, as discussed above. The D ^ D correlation function with a nonzero ffnal 
state momentum is used to extract w via Eqs. (2.6) and (2.7). 
The analysis below also requires the two-point function 

C^{t,p) = E e'"" (^xit, ^)Ox{0, 0)), (3.4) 

X 

where X = B ot D. 

We simulate the daughter meson with two different choices for the interpolation operator 
Ox'- with a IS-wave smearing and without any smearing [42]. The smearing optimizes 
the overlap of the operator with the ground-state wave function of the meson, so the two 
choices have different excited-state contributions but the same energies. For the three-point 
correlation functions, we always use a IS-smearing source for the extended quark propagator. 

B. From correlators to form factors 

In general, two- and three-point functions, such as those in Eqs. (3.1)-(3.4), can be 
expressed as [43], 

C^{t,p) = ^(-l)^-*|Z,(p)|2 [e-^^(p)* + e-^^^(P)(^-*)] , (3.5) 

k=0 

C3^^"^(0,t,T;p) = ^^(-l)^-*(-l)^(^-*)A^^(p) e-^^(P)* e-*^^(^-*) (3.6) 

k £ 

where Ek [Mg) are the energy levels of Y (X) and A'^^^ are coefficients of the transition 
Xe — Ifc. If the time differences between the source (0) and vector current (t) and that 
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between vector current and sink (T) are sufficiently large, i.e., \t\ — > oo and \T — t\ 
only the lowest energy level will survive. Then, we have {Bg — > Dg follows similarly) 



oo, 



hi ^ 



C3^^*^(0, t, T; 0) C3^^'^(0, t, T; 0) ' 



(jDVW 



(0,t,T;p) / |Zo(0)| Eoip) 



C3^^'^(0,t,r;p) 
Cr'^{0,t,T;p) 

C.^y'^{0,t,T;p) 



^[Eo{p)-Eo{0)]t 



|Zo(p)| V ^o(O) 



(3.7) 

(3.8) 

(3.9) 
(3.10) 



where means that the left-hand side is output of an analysis procedure. In practice, the 
separations between the current insertion and the source/sink, t and T — t, are often not 
large enough to suppress the excited states completely. The factor inside the parentheses 
in Eq. (3.8) cancels the leading time dependence of the ratio of three-point functions with 
different momenta; both Zq{p) and Eq{p) come from fitting C2 as suggested by Eq. (3.5). 
Instead of fitting for plateaus, we extract the matrix-element ratios on the left-hand sides of 
Eqs. (3.7)-(3.10) by fitting the right-hand sides in a way that incorporates excited states. 
We can write the ffist few terms in Eq.(3.6) as 



^02 (^) + ^20 (^) + higher excitations. 



(3.11) 



where A'j^^{t) = A'^^ Q-Ekt^-MeiT-t) _ rpj^^ terms Aiq, Aqi, An are the contributions from the 
opposite-parity states that are introduced by the operator involving staggered quarks. We 
can reduce their effects by making use of the fact that they oscillate with either t or T. 
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FIG. 2. A sample fit of the single ratio di with the final D meson carrying a momentum p = 
27r(l, 0, 0)/L on the ensemble C020. The left graph shows the data with IS-smearing (filled circles, 
data along the bottom curve) and that without smearing (open circles, data along the top curve) 
fitted separately (solid curves) and simultaneously (dashed curves). Fit results over the interval 
t G [^miii) 10] are shown in the right graph and are seen to be stable for tmin > 1. 
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The contamination from Aiq and Aqi is minor because no obvious oscillation is visible with 
any of the three-point functions. Based on the plots of the ratios calculated from different 
source-sink separations T, there is no sizable contribution from An either. With the four 
source times, which are evenly distributed in the temporal interval of the lattice, we separate 
two of the sinks from the sources by an even number of time slices and the other two by a 
neighboring odd number. We apply the following averaging method to further reduce the 
effect from An [42] 



1 



1 



1 



i?(0, t, T) = -R{0, t, T) + -i?(0, t, T + 1) + -i?(0, t + 1, T + 1) 



(3.12) 



where R stands for any of the correlation-function ratios corresponding to R^, ai, bi, or c?,. 
With this averaging procedure, the An terms are suppressed by a factor of 6-10 and the Aqi 
and Aio terms, already small, are suppressed by a factor of about 2. Hence, the systematic 
error arising from neglecting terms Aqi, Aiq and An can be safely dropped. 

The foregoing analysis enables us to use a simple fitting scheme for the ratios including 
only the contributions from Aqq, A20, Aq2. At the lowest order, the functional forms for 
determining Oj, bi and di are then 



Cr*''{0,t,T;O) \Zo{p 
= di [l 

bi [I 




(0,t,T;p) 



Ci>^'^{0,t,T-p 



(JDV'D 



[I 



02 e 



9o2 e 



^02 e 



AM{T~t) 



+ 



/20 e-^^(^)* + 



20 e 



-A£;(0)tl 5t 



-AM(T-t) _|_ 



~AE{0){T~t) 



ho e 



~AEip)tl 



^20 e 



~AE{p)t 



(3.13) 
(3.14) 

(3.15) 



The fit parameter 5 in Eq. (3.13) allows for imperfect cancellation of the leading t depen- 
dence. The parameters AE{p) = E2{p) — Eo(p) (AM = M2 — Mq) denote the splittings 
between the D-meson energy (i?-meson mass) and its first radial excitation. We find that 
the double ratio R^ is so weakly affected by excited states that it suffices to fit it to a 
constant in t. Adding terms to describe excited states in R_^. leads to changes no bigger 
than the statistical errors. The energy splittings AE{p), AM in these expressions are also 
constrained by the two-point functions, Eq. (3.5). 

We employ two fit procedures to determine the ratios. One, which we explain first, is 
simpler. The other is more complicated but yields better results, as we explain below, so we 
take it as our primary analysis and use the simple method as a cross check. 



C. Two-step Fit 

The simpler method proceeds in two steps. We first fit the two-point functions to obtain 
the energies, energy splittings, and overlaps: E{p), AE{p), AM, and Z{p). We use con- 
strained curve fitting and priors [44, 45]. The priors in these fits impose no real constraint. 
Second, we take the energy splittings from these two-point fits as priors when fitting the 
ratios of three-point functions. The priors for the amplitudes on the right-hand sides of 
Eqs. (3.13)-(3.15) are again taken wide enough to impose no real constraint. Below we call 
this approach the "two-step fit." 
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As discussed above, we have data for smeared and local interpolating operators. These dif- 
ferent correlator ratios have different excited-state contributions. To determine the matrix- 
element ratios, we fit the correlation-function ratios of the two smearing types either sep- 
arately or simultaneously, as illustrated in Fig. 2. Although the two smearing types follow 
rather different curves, they arrive at consistent values of Oj, bi, and di as described in 
Eq. (3.13)-(3.15). Figure 2 shows that the separate fit and the simultaneous fit give consis- 
tent results for the ratio d^. The other two single ratios Oj, 6, can be determined in a similar 
way, which can be seen in the sample fit in Fig. 3. 



D. Combined Fit 



Our preferred fit treats the two-point functions and three-point-function ratios simulta- 
neously. This approach allows all correlation functions to influence the output energies and 
overlaps, combining all information at hand, including the correlations between the two- 
and various three-point functions. We refer to this procedure as our "combined fit" , in con- 
trast with the two-step fit described above. In particular, all correlation functions are then 
treated on the same footing in the determination of the energy splittings. Figure 4 shows 
the relationships between the ratios and two-point functions, building up constraints among 
correlation functions of zero and nonzero final momenta. The single ratios also require 
the two-point Z{p) factors from the relevant three two-point functions. 

The combined fit is repeated for each value of the momentum. The resulting single 
ratios a^, bi, and di and the double ratio R+ at the corresponding recoil are then used to 
compute the form-factor ratios h±{w)/h^{l) and To verify the results from the two 

approaches, we compare the resulting h+{w) of the two coarse ensembles C007, C020 in 
Fig. 5. In general, they are in very good agreement at the five values of w where we have 
data, while the combined fit has slightly better precision at small recoil. 

The combined fit turns out to have other more important advantages over the two-step 
fit. First, the resulting hj-{w) / h^{l) is more stable with the combined fit than the two-step 
fit, where the fitting range must be determined carefully. This stability stems from the fact 
that the combined fit does a better job of resolving the correlated statistical fluctuations at 
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FIG. 3. A sample fit of the three single ratios Oj, bi, and di from the ensemble C020. The circles 
and triangles correspond to the data without smearing and with IS" smearing, respectively. The 
final D meson carries a spatial momentum p = 27r(l, 0, 0) /L. The dashed curves indicate the best 
fits and they are in good agreement with the data points. The horizontal lines show the resulting 
ratios ai, bi, and di, as well as the ranges included in the fits. 
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zero and nonzero recoil. Second, the combined fit helps to reduce the systematic error due to 
excited states. We account for the excited state contribution in the fit for single ratios to work 
around the fact that sink-source separations cannot be taken satisfactorily large. Although 
the towers of the excited states are not the focus of this paper, the proper accounting of 
their contributions is important, because they influence ratio results. The combined fit 
procedure resolves the excited states more stably than the two-step fit. Third, the resulting 
form factors h±{w) using the combined fit at these values of w are more consistent with 
each other than those using the two-step fit. This can be seen when one attempts to fit the 
results at different recoil to a chiral effective theory. Section IV shows that the combined 
fit procedure results in a more reliable chiral extrapolation than the two-step fit, and hence 
we use it throughout our analysis. That said, as shown in Fig. 5, the two fitting procedures 
give seemingly good pointwise consistency. 

Although the combined fit method helps in many aspects, it requires the handling of a 
larger data set and correspondingly a much larger covariance matrix. We use the jackknife 
method with single elimination to calculate the covariance matrices because the data samples 
show a very small autocorrelation time (less than 1). To reduce the time searching for the 
minimum of x^? we take the output of the two-step fit as the initial guess for the combined 
fit. It is worth mentioning a small complication. When fitting the single ratio a^, we need in 
advance both the ground-state energy Eq and the wave function renormalization factor Z to 
suppress the time dependence of the ratio. With the combined fit procedure, Eq and Z are 
refined through the two-point functions which are part of the combined fitting. So in the 
actual analysis, we take the results of Eq and Z from the combined fit and plug them back 
to suppress the time dependence of the three-point function for Oj. We need to iterate such 
a process a few times until the fitting results stabilize. We find that this iteration converges 
within two or three steps. 




FIG. 4. Diagram showing which correlators and correlator ratios influence energies and amplitude 
ratios. Energy splittings of the initial and final mesons B(0), D{0), and D(p) are determined from 
the two-point functions (boxes) as well as the ratios (circles). The lines connecting them indicate 
their common dependence on the splittings. Altogether seven correlation functions are included in 
the combined fit. 
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IV. CHIRAL- CONTINUUM EXTRAPOLATION 



Given the light-quark masses in Table I, we extrapolate the results to the physical value 
guided by chiral effective theory. In the case of h+{w), for a heavy-light system, we follow 
rooted staggered chiral perturbation theory (rS^PT) [46-50]. The specific application to 
the case of 5 — D'-*^ at zero recoil is provided in Refs. [42, 51]. The continuum xPT 
for the semileptonic B — )■ D^*^ form factor at nonzero recoil has been derived at next-to- 
leading order (NLO) in Ref. [52], and the generalization to rS^PT for B ^ D is given 
in the Appendix. In the case of h^{w), the leading correction is simply a constant that 
is inversely proportional to the charm quark mass. To describe the simulated data, we 
also must parametrize the recoil dependence to quadratic order around zero recoil. The 
expansion coefficients are related to the slope and curvature of the form factors. 

Thus, we can write the general expression for h±{w) with NLO rS^PT and higher-order 
analytic terms incorporating lattice-spacing dependence as 

hf{w) = 1 - pKw - 1) + Kiw - ir + ^±i^ + ^^logSi,„„p(A„^) 

c J 

+ Co,+ + ci^+ {2mi + rrih) + Ca,+a^ (4.1) 

h^^\w) = — - p'i{w - 1) + k_{w - 1)^ + Co - + ci _ (2m; + rrih) + Ca .a^ (4.2) 
rric , , , 

where — and 2k± are the slopes and curvatures of the form factors, while X± are low- 
energy constants. In the case of h+, depends on the chiral scale A,^ in such a way as to 
cancel the A^ dependence of the nonanalytic terms ("chiral logs"). These terms, denoted 
here as logs^.j^op, are given by the terms appearing inside the square brackets in Eqs. (A. 2) 
and (A. 3) of the Appendix. The other form factor h_{w) has no nonanalyticity at one loop. 
Equations (4.1) and (4.2) also include terms depending linearly on the valence {rrix) and sea 
{mi and m^) quark masses, with coefficients Co(i),±. These terms are next-to-next-to-leading 
order in the chiral expansion and are needed to describe the data with rrix or mi ^ ^mg. 
Generic lattice-spacing dependence is described by the terms with coefficients Ca^±- 
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FIG. 5. Comparison of h^{w) obtained from the two-step fit and combined fit for the two coarse 
ensembles C020 (left) and C007 (right). 
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TABLE II. Input parameters for the chiral extrapolation [22, 53]. 
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0.2052872 
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rfa^AT 


0.3268607 


0.3268607 


0.1153820 


0.1153820 


rfa^Av 


0.4391099 


0.4391099 


0.1523710 


0.1523710 


rfa^Aj 


0.5369975 


0.5369975 


0.2062070 


0.2062070 




-0.05 


-0.05 


-0.03 


-0.03 




-0.28 


-0.30 


-0.15 


-0.16 



We treat the chiral extrapolations of the B ^ D and Bg — t- data slightly differently. 
For B ^ D, we analyze only full QCD data points, i.e., = mi. Then, since the strange 
sea quark in all ensembles is tuned within several per cent of its physical mass, ruh ~ m^, the 
dependence of the form factors on the sea and valence quark masses cannot be disentangled. 
Therefore, we drop the parameter ci^+ when fitting the B ^ D data. For our Bg — )■ Dg 
data, on the other hand, the strange valence quark is tuned close to its physical mass for all 
the ensembles we analyze, rrix = mg. As a result, we cannot disentangle the valence quark 
dependence. Therefore, we discard the parameter co,± when fitting the Bg — > Dg data, and 
estimate the tuning error of rrig a posteriori in Sec. VI D. 

The rSxPT expression for logS]^_ioop(A^, w) contains several low-energy constants used in 
rSxPT to describe the masses and decay constant of light pseudoscalar mesons. The values 
we use for these parameters are taken from Refs. [22, 53] and given in Table II. The xPT 
expressions also require the D(^gyD'^^-^ splitting A^'^^ and the pion decay constant we take 



B->D: x7dof=5/20,p-value=1.0 



B^->D^: 5c7dof=7.6/20, p-value=0.99 





FIG. 6. Chiral-continuum extrapolation of h^{w) for B ^ D (left) and Bg — )• Dg (right) decays 
based on the four ensembles listed in Table I. The blue bands show only the statistical errors and 
the red curves are the chiral and continuum limits. 
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TABLE III. Chiral extrapolation of the two-step and combined fit results on the coarse ensembles. 

Form Two-step fit Combined Fit 

factor x^/d-o.f. p value x^/d-o.f. p value 

hf-^^{w) TT/T 014 3.7/7 0.81 

h^s->D''(w) 11/7 0.13 5/7 0.68 



both from Ref. [15]. To combine data from both lattice spacings, we convert dimensionful 
quantities to ri units, where ri is the distance defined via the interquark force by r^F(ri) = 1 
[54, 55]. We take n/a from Refs. [22, 53]. 

Unfortunately, the D*-D-'k coupling gD*D-K and similar couplings with strange mesons, 
which appear in the coefficient gjj*j^^/167r'^f^ of the chiral log terms, are not known with 
good precision. We appeal to various estimates of gD*Dn available in the literature, in- 
cluding CLEO's measurement of the D* width: gD*Dn = 0.59(7) [56]; quenched lattice 
QCD: gD*D-K = 0-67(8) (Ig) [57]; a fit to various experimental data, including the D* width: 
go'D-K = 0.51 (no error reported) [58]; two-flavor lattice QCD in the static limit: gD*Dn = 
0.516(51) [59]; and 2+1-flavor lattice QCD in the static limit: gD-'Drr = 0.449(51) [60]. In this 
calculation, we include gD*Dn as a parameter in the constrained fit with a prior 0.51 ± 0.20. 

In Sec. Ill, we compared the two fitting procedures, two-step fit and combined fit, with 
which we obtain the single ratios a^, 6j, di and the double ratio R^. At each w where we have 
data, h^{w) from the two procedures are in good agreement (within la). We then fit the 
resulting h+{w) for the coarse ensembles (C020 and C007) to Eq. (4.1) without the analytic 
terms and the dependence (NLO). The results are shown in Table III. It is apparent 
that the results from the combined fit procedure are better described by the chiral effective 
theory that we employ, giving a x^/d-o.f. = 0.53 for B ^ D (compared to 1.6 from the 
two-step fit). A similar observation can be found in the case of Bg — )■ Ds- This indicates 
that the correlations among the ratios and those among different kinematic points are better 
resolved by the combined fit, and hence we follow this procedure for the entire analysis. 



B->D: x^/dof=8.3/ll,p-value=0.69 Bj,->Dj,: 5C^/dof=17.6/l 1, p-value=0.09 




FIG. 7. Chiral-continuum extrapolation of h-{w) for B ^ D (left) and Bg — )• Dg (right) decays 
based on the four ensembles listed in Table I. The blue bands show only the statistical errors and 
the red curves are the chiral and continuum limits. 
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The results of the chiral- continuum extrapolation of hj^{w) ioT B ^ D and Bg — t- Dg are 
plotted in Fig. 6. With the large number of configurations we have for the four ensembles, we 
are able to determine the form factors hj^{w) with statistical errors at the level of ~ 0.5% at 
zero recoil, increasing to ~ 1.5% at w = 1.15. The form factor h+{w) for both of the B ^ D 
and Bg — )■ Dg decays exhibits a small dependence on the light-quark masses and lattice 
spacings, so the extrapolated physical values are close to the lattice data. The difference 



between h^~^^{w) and h^"^^''(w) is also small. The form factor with strange spectator 
h^''^^"{w) shows a steeper slope and larger curvature — = 1.26(09) and /c+ = 1.15(9) — 
than its 5 ^ D counterpart— = 1.14(10) and k+ = 0.87(13). 

The results of the chiral-continuum extrapolation of h_{w) for B ^ D and Bg — t- Dg 
are plotted in Fig. 7. Here light-quark mass (both sea and spectator) and lattice spacing 
dependence are visible. We find an k, 0.04 difference between the two /i_ at zero recoil. From 
Eq. (2.11), however, this effect does not cause much difference in /q. We therefore anticipate 
that the ?7-spin symmetry breaking effect is smaller than what was found in Ref. [21]. Such 
an observation is in accord with the recent lattice calculations of /+, /o in the -D(<,) — )■ vr(ii') 
decays [25]. 

With the results from the chiral-continuum extrapolation in hand, we now convert h±{w) 
into f+{q^) and fo{q^) with Eqs. (2.10) and (2.11) and the physical B and D masses [15]. 
To gain confidence in our procedures, let us compare the resulting /+ with experimental 
measurements. The differential decay rate oi B ^ D is given by 



MI{Mb + Mof{w^ - lf'^\V,,\''\G{w){' 



dw 487r3 
where is the Fermi constant, and it is conventional to introduce 



(4.3) 



1 + r 



(4.4) 
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FIG. 8. The slope of the form factor G{w) from the chiral extrapolation is compared with var- 
ious experimental measurements from Belle [61], CLEO [62], BaBar (tagged) [63], and BaBar 
(global) [64] respectively. 
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Experiments report the zero- recoil form factor and the relative form factor slope 

= —Q'{1)/Q{1) [65]. From our extrapolated data, we find ^(1) = 1.058(9)stat., which is 
consistent with the previous unquenched lattice-QCD result 1.074(18)stat.(16)syst. [66]. The 
measured slope is related to parameters of our chiral extrapolation via 



1 



P 



r + 1 



P- 



(4.5) 



where pf^gg is the slope of the NLO logarithm at zero recoil. In Fig. 8, we compare the slope 
of the B D form factor at = 1 with experiment. We find = 1.25(5)stat.5 where the 
error is obtained from 300 bootstrap samples. This value is in good agreement with the 
experimental results from Belle, CLEO and BaBar [65] and the average of these from the 
Heavy Flavor Averaging Group, = 1.18(6) [67]. 

Note that a determination of Q{w) with full error budget is beyond the scope of this 
paper. A comprehensive effort to do so is in progress [68]. Here we are merely satisfied to 
see that the main ingredients of our analysis are compatible with the available experimental 
data. 



V. z PARAMETRIZATION 

To minimize discretization effects, the final D meson momentum p should not be taken 
too large, so the calculations are restricted to small recoil, w < 1.17. However, the form- 
factor ratio that we are trying to compute ultimately needs to be evaluated near maximum 
recoil, w ~ 1.6, which appears to require a considerable extrapolation. Fortunately, the 
extrapolation can be guided by the model-independent z parametrization [23, 69, 70]. As 
shown, for example, in Ref. [71], this strategy is effective for extrapolating lattice-QCD data. 

One introduces the variable 

ziw) - ^^^"^ (5 1) 

which maps the physical domain of the form factors into the unit disk. The form factors 
can be expressed as 

^ oo 

where P{z) and (j){z) are called, respectively, the Blaschke factor and the outer function. 
The range of w for B^D is l<w< 1.589, which corresponds to < z < 0.0644. With the 
choice of outer functions given below, unitarity sets a bound on the expansion coefficients. 



„r < 1^ (5.3) 



Because of this constraint and the restricted range of z, the expansion in Eq. (5.2) converges, 
and one can parametrize the form factors with only a few terms. In this analysis, we truncate 
the expansion at the 2;^ term, which is enough for our data. 

Although we are primarily interested in /q, we apply the z expansion to /+ and /o 
simultaneously, incorporating the kinematic constraint /o(0) = /+(0). We take the outer 
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functions to be 



^o(^) = $0 (1 + z) (1 - zf'^ [(1 +r){l-z) + 2v^(l + z)] ' 
^.+ (^) = $+ (1 + - [(1 + r){l-z) + 2v^(l + z 



(5.4) 
(5.5) 

where we choose the constants $o = 0.5299 and = 1.1213 such that it matches the same 
normahzation as in Ref. [69]. These exact values can be combined with the fit resuhs, given 
below, to reconstitute the form factors. 

The purpose of the Blaschke factor is to remove poles outside the physical region, i.e., 
ai w <1. In general, 

z{w) - z{wp) 



P z) 



nr 

p=i 



- z{w) z{wp) ' 



(5.6) 



where z{wp), with Wp = [1 + — M^/M]^ ^]/2?", marks the pole position, and the product 
can run over states with Mp < Mb^^-^ + ^d^^) ■ In the case at hand, these poles appear at 
the masses of = 1~ vector {J^ = scalar) mesons for /_|_ (/q). 

We need to put the results of the chiral-continuum extrapolation in the 2;-parametrized 
form of Eq. (5.2). We do so by generating synthetic data from the chiral-continuum ex- 
trapolated curves and fitting for the corresponding parameters in Eq. (5.2). Our chiral- 
continuum fit has eleven free parameters, including ones to describe the dependence 
in h±{w). Such terms vanish in the continuum limit, leaving only nine physical parameters. 
Furthermore, due to the small contribution from the terms of go* d-k and Co(i),+ and due to the 
correlation between h±{w), we end up with effectively six free modes in the synthetic data. 



>D z-expansion: % /dof=0.31/l, p-value=0.58 



I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 




B5->Dg z-expansion: % /dof=0.96/l, p-value=0.33 
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FIG. 9. z expansion of form factors to the maximum recoil point for both B ^ D and Bs Dg. 
The diamonds and circles for z < 0.02 are the synthetic data derived from the chiral-continuum 
extrapolation, and the solid curves and error band are the results of the z expansion. The dashed 
curves show how the chiral-continuum extrapolation extends into the region where the extrapolation 
may not be less trustworthy. The points near z = 0.06 correspond to the desired = and 
M|.; squares (triangles) correspond to z fits with (without) the B* pole in the Blaschke factor. 
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Thus, we generate the synthetic data at three evenly-spaced values w G {1.0, 1.08, 1.16} in 
the region where we have data. In Sec. VI B, we show that the final form-factor ratio is not 
sensitive to the details of these synthetic data. We start with the trivial Blaschke factors (no 
poles) and fit the B ^ D and Bg — )■ Dg form factors up to the term in Eq. (5.2). The re- 
sults of these fits are shown in Fig. 9. To examine the effect of poles on the shape of the form 
factor /+, we also try a fit with a one-pole Blaschke factor at mass Mp = Mb* = 6.330 GeV, 
where this lowest B* is a prediction of lattice QCD [73]. This Blaschke factor affects the 
extrapolated results with a deviation of about 0.3%. The dashed lines in Fig. 9 show the 
extrapolation based solely on the chiral fit, showing that the z expansion plays an important 
role in controlling the total error. 

We perform the ^-expansion fit without constraints on ag, Oi, and a2, and setting the rest 
to zero. In Sec. VI B, we discuss fits with more parameters, constrained then by Eq. (5.3). 
We constrain the fit with the relation /o(0) = /+(0), by demanding |/o(0) — /+(0)| < 6 where 
6 can be chosen arbitrarily small. Once 6 is small enough, its actual value has no effect on 
the fit. 

To check the form factor shape obtained from the z expansion, we can compare the 
B ^ D decay with the latest published measurement from BaBar [72]. The comparison is 
shown in Fig. 10, with \Vcb\ = 41.4 x 10~^ as used in Ref. [72]. As one can see, the shape 
of our form factor agrees well with experiment over the full kinematic range. As above, we 
note that this comparison is made without a full treatment of the systematic errors on Q{w). 
A thorough treatment with full error analysis, aimed at determining \Vcb\, will be covered 
elsewhere; see Ref. [68] for a progress report. 

For completeness, we give the results of the z fit in Table IV. The correlation matrix does 
not include full systematics, but with this information, the reader can reproduce the curves 
and error bands in Figs. 9 and 10. Note, however, that the parameters of the nonstrange 
and strange form factors are also correlated. 
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FIG. 10. Comparison of the form factor shape of ^(w)|V^fe| with BaBar's measurements [72]. 
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TABLE IV. Best-fit values and correlation matrix p^n of the simultaneous 3-term z expansion 
of /+ and /o, with statistical (post extrapolation) errors only. Top: B ^ D; bottom: Bg — )■ Ds- 
Note that the fit parameters are correlated between the B ^ D and Bg — )■ Dg processes. 
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VI. SYSTEMATIC ERRORS 

We now discuss the systematic errors in our analysis. Owing to the similarity between 
the B ^ D and Bg — )■ Dg processes, the systematic errors in the ratio of the form factors 
largely cancel, by design. To assess the systematic uncertainties, we have repeated the chiral- 
continuum and z extrapolations with different choices. The values of fQ^\M^), flf\M]^), 
and flf\M^) / flf'\M'^) resulting from these variations are listed in Table V. We summarize 
the final error budget in Table VI. As our standard analysis, we use Eqs. (4.1) and (4.2) 
for the chiral-continuum fit, dropping ci^± ioi B ^ D and co,± for Bg — )■ Dg. We fit the 
coupling gD*DiT using a constrained fit [44] with the prior 0.51(20). We take the synthetic 
data points a.t w = 1.0, 1.08 and 1.16 and include the trivial Blaschke factor (no poles). We 
use ri = 0.3117 fm to convert the necessary physical inputs (like /^r) to ri units. 

A. Chiral extrapolation 

Our chiral extrapolation is based on the rS^PT formalism shown in Eqs. (4.1) and (4.2) 
and the Appendix. The systematic errors arising here can be divided into two categories: 
the error associated with the one-loop contribution itself and the error associated with the 
partial inclusion of NNLO analytic terms. Throughout our analysis, we keep the slope 
and curvature k, because they determine the w dependence of the form factors. 

An uncertainty in the contribution from the NLO logarithm stems from the uncertainty 
of the D*-D-T[ coupling gD*D-n- Although the NLO logarithms to hj^{w) make a small 
contribution when evaluated with the quark masses for which we have data, of order 10~^ 
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TABLE V. Values of the form factors /q ' and their ratio for several variants of the fitting 
procedure. The second panel shows the results for different choices of chiral extrapolation fit 
function. Note that for the B ^ D form factor we use the fit functions labeled "val" while for the 
Bs —5- Ds form factor we use the fit functions labeled "sea" in Eqs. (6.1)-(6.5). The third panel 
shows the results for different choices of the z expansion. The final panel shows the results for 
different values of parametric inputs: the lattice scale and light- and strange-quark masses. 
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i(;-Range[l,1.20] 


0.646(19) 


0.618(33) 


1.048(45) 


i(;-Range[l,1.24] 


0.653(19) 


0.623(34) 


1.049(46) 


Truncated at 


0.632(22) 


0.607(36) 


1.042(45) 


Truncated at 


0.632(22) 


0.608(36) 


1.043(46) 


ri = 0.321 fm 


0.638(18) 


0.611(32) 


1.047(44) 


rus Ic shift 


0.638(18) 


0.611(32) 


1.046(44) 


mi la shift 


0.639(18) 


0.612(32) 


1.046(44) 



at our lightest simulated quark mass, the logarithm affects the form factor shape of h^{w) 
through its w dependence. Thus, the uncertainty in Qd'Otx becomes more important as we 
extend our calculations to large recoil. We include Qd'O-k in the chiral-continuum fit with 
the prior 0.51(20), which describes the data well, so we do not assign an additional error 
due to the uncertainty of Qo'D-k- 

We now look at variations from fitting with and without the NNLO analytic terms. For 
ease of discussion, let us break the chiral fitting scheme into the following different pieces: 



NLO+ = 1 - pI{w - 1) + k^{w -lf + -± + J|!i^logs,.i„,p(A^,^), (6.1) 

^LO- = l-p'_{w-l) + k_{w - 1)2 + — , (6.2) 

pNNLO± = NLO^ + Co,±m,, (6.3) 
pNNLO± = NLO± + ci,±(2m^ + m,), (6.4) 
pNNLOi„,_/,,, = pNNLOi_/,,, + c.,±a^ (6.5) 

where "pNNLO" stands for partial NNLO, because we include only analytic terms. We are 
not aware of any full NNLO calculations with nonzero final D-meson momentum. 

As mentioned in Sec. IV, the form factor /i+(w) shows a weak dependence on the light 
sea and valence quark masses, and the data are already well-described by NLO+, with 
xVd.o.f. = 0.68 and 0.83, respectively, for h^-^^{w) and hl'^^'{w). Adding the pNNLO+ 
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TABLE VI. The error budget of the form-factor ratio discussed in text. 



Source of error 




Statistics © chiral extrapolation 


4.2% 


z expansion 


0.6% 


Scale r\ 


0.1% 


Mistuned nis 


0.1% 


Mistuned m/ 


0.1% 


Heavy-quark [k) tuning 


0.6% 


Heavy-quark discretization 


1.0% 



terms, the hj^ fits remain good. As seen in Figs. 6 and 7, the data for h-{w) exhibit a more 
significant dependence on the lattice spacing and on the sea- and spectator-quark masses. 
Unsurprisingly, the NLO~ fit of h^{w) leads to large x^/d-o.f. = 1.7 and 1.9, respectively, for 
j^B-^D g^^^ j^Bs^Ds^ rjij^g NLO xPT correction to h^{w), denoted X-frric, does not depend 
on quark mass or lattice spacing, so the observed dependence in the data must be described 
by pNNLO terms. Moving through the pNNLO functional forms in Eqs. (6.3)-(6.5), we 
find that the pNNLO;^^^^i fit of h^--*^{w) improves nicely, xV^-o.f. = 0.93, but pNNLO" 

fit of h^''^^''{w) less so, x^/d-o.f. = 1.8. These can be contrasted with the standard fits, 
pNNLO;!;^^^! for 5 -> L) and pNNLO^^,^^ for D„ with h_ xVd.o.f. = 0.75 and 1.6, 

respectively. As seen in Table V, these less good fits all lie well within the extrapolated 
statistical error of the standard fits. We therefore treat these alternatives as cross checks 
and do not add an additional error here. 



B. z expansion 

Although the z expansion provides a model-independent parametrization of the form 
factors /o and /+, the final results may depend on three kinds of choices made within 
this framework. First, the expansion coefficients may depend on the number and range of 
synthetic data points. Second, the shape of the form factor may be affected by the number 
of poles in the Blaschke factor, particularly for f^. Last, the shape may be affected by the 
truncation of the series in z if one does not include enough terms. 

To estimate the uncertainty from the synthetic data, we vary the w range over which they 
are generated. We repeat the z expansion with w in the intervals [1, 1.08], [1, 1.12], [1, 1.20], 
[1,1.24], choosing three evenly-spaced points for both form factors. We also try fits with 
more synthetic data than underlying parameters, in which case some of the information is 
spurious, leading to tiny eigenvalues in the synthetic-data covariance matrix. We remove 
the corresponding mode(s) with singular-value decomposition. Although the form factors 

flf''^\o) each vary with these alternative choices by about la, the ratio flf\M^)/flf'\M'^) 
is negligibly affected; cf. Fig. 11. We take the maximum deviation from 1.046, which is 
0.004, as the systematic error. 

To estimate the uncertainty from the poles in the unphysical region, we repeat the z ex- 
pansion fit by including a B* pole in the Blaschke factor for /+. We take Mb* = 6.330 GeV 
from a lattice-QCD calculation on the MILC ensembles [73]. Recall that /+ infiuences /o 
near maximum recoil via the kinematic constraint /o(0) = /+(0). We find that the effect is 
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rather small, as shown in Fig. 11, leading to a difference of only 0.3% in the ratio of /q. 

To estimate the uncertainty from the truncation of higher order terms in the z expansion, 
we perform the z expansion fit by including and terms and incorporating the unitarity 
constraints on the coefficients, z.e., < 1. As can be seen in Table V, both the form factors 
and their ratios stabilize when higher order term z^ (or further z^^ is included. This results 
in a 0.3% difference in the form-factor ratio. 

The total systematic error including all these effects added in quadrature is 0.6%. 



C. Scale ri dependence 

As discussed in Sec. IV, we convert our data to r\ units with ri/a from Ref. [22]. To 
convert to physical units, we must choose a value of in physical units. Our choice is 
based on MILC's analysis of /^r, which leads to ri = 0.3117(6) (;1^3^) fm [22]. To estimate 
the error, we also consider an early value from HPQCD based on the 2S'-1S' splitting of the 
T resonances, = 0.321(5) fm [74], and repeat the whole analysis with this value. We 
find a negligible shift, ~ 0.1%, in the form-factor ratio, because of the cancellation between 
Bs Ds and B ^ D. 



D. Light- and strange-quark mass dependence 

The physical light and strange quark masses are determined from the analysis of light 
pseudoscalar meson masses and decay constants [22, 75]. To estimate the error, we repeat 
the chiral-continuum extrapolation varying the masses by ila. We do so twice, once varying 
the physical light quark mass by 3.1% for B ^ D; and again varying the physical strange- 
quark mass by 3.4% for Bs — > Dg. In both cases, we find a shift on fQ\M^)/flf'\M'^) less 
than 0.1%, which is much smaller than other errors in this analysis. 



1.15 



"1 \ 1 1 r 



o 



1.1 - 



1.05 - 



1 - 



O O C) 



0.95 



J \ I I L 



1.08 1.12 1.16 1.20 1.24 f+-pole 



FIG. 11. Systematic error due to the z expansion. The open circles show the effect of varying the 
synthetic data used, and the filled square shows the effect of adding a pole to the Blaschke factor. 
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E. Heavy-quark mass dependence 

We have not generated data for a wide-enough range of Kb and Kc to determine directly 
the heavy-quark mass dependence of the Bg — ?■ Dg form factors. Since our main focus is a U- 
spin breaking ratio, we rely instead on other such ratios computed on the same ensembles, in 
particular the form factor for B — )■ D* at zero recoil /iAi(l) [42, 76] and the ratio of leptonic 
decay constants fss/ fB+ [77]. 

In the case of which is very similar to we find the K-tuning error to be 

0.56% of /lAi(l) and 4.8% of 1 - hA,{l) [76]. In the case of ^/ = fBjfB+, we find the 
K-tuning error to be 0.41% of ^/ and 2.2% of .^j — 1 [77]. The first of these four estimates 
yields the largest absolute error on ft^/fl^'^\ namely 0.6%. This error estimate is still much 
smaller than the overall error in this analysis. 

F. Heavy-quark mass discretization and matching 

We also use our work on /iAi(l) [42, 76] to guide and estimate heavy-quark discretization 
errors, both power-law and radiative effects. For h^-^^l), we find a 1.0% error from discretiza- 
tion effects, and a 0.3% error from matching. Since the present calculation matches only at 
tree level, the corresponding errors here are order instead of a^. For a ?7-spin-breaking 
ratio such as ours, the discretization error is further suppressed by (m^ — md) / Aqcd- Since 
{nis — niii) /{as Aqcd) ~ |, there is not much change. From the structure of Eq. (2.11), the 
matching error stemming from is asiml — m^)a^, which is negligible, but the matching 
error from h_ leads to an error on /q of order as{h_/h+){ms — tridj/Kq^cxi ~ 0.5%. Taking 
a 1% error for these effects seems reasonable yet does not influence the total error budget 
much. 

G. Summary 

Let us now summarize our results. Table V hsts the values of f^"\M^), f^'^\Ml) and 
their ratio flf\M^)/f^'^\M'^) under the variations in the analysis explained above. The 
resulting error budget is given in Table VI, based on which, we arrive at our final result 
given in Eq. (1.6). The systematic error is the sum of the listed systematic errors added 
in quadrature. Shifting the argument of the denominator slightly and following the same 
analysis steps, we obtain Eq. (1.7). 

VH. CONCLUSION AND DISCUSSION 

To conclude, we provide the first lattice-QCD calculation of the form-factor ratio 
f^'\M^)/f^'^\M^). Our result leads to the factor 7V> = 1.094(88)(30), which is signif- 
icantly closer to unity than the sum-rule estimate [21], Afp^ = 1.24(8) (or Afp^ = 1.3(1) 
[17]) used in previous analyses of hadronic fs/ fd [18, 20]. As noted above, the lack of 
significant ?7-spin breaking observed in this calculation is in accord with other lattice-QCD 
calculations of similar form factors [25]. 

We now examine how our new value of JVp affects the fragmentation-fraction ratio fs/ fd- 
LHC6 measures fs/fd via BR(i?° — D+7r")/BR(i?° — D~^K~), using the sum- rule estimate 
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J\fp^, and finds fs/ fd = 0.250(24)stat(17)syst(17)thco [20]. Since Mp is not correlated with any 
other quantity in Eq. (1.3), we easily find that the fragmentation ratio should become 

^ = 0.283(27)stat(19)syst(24)theo, (7.1) 
Id 

where the errors have also been scaled accordingly. Superficially, our theoretical error is 
slightly larger than that obtained with the sum-rule estimate — 8.5% vs. 6.5%. Our error, 
however, is straightforward to improve, since it is dominated by Monte Carlo statistics, 
propagated through the chiral-continuum and z extrapolations, as seen in Table V. 

Fleischer, Serra, and Tuning have proposed a second hadronic approach based on the ratio 
BR(5° L)+7r")/BR(i?° D+vr") [18]. A comphcation is that a l^-exchange diagram 
also contributes to the — > D^ir" decay, leading to an additional factor N'e in the analog 
of Eq. (1.3). It is estimated to be Me = 0.966(75) [20]. This method requires a similar input 
of the form-factor ratio A^^ = [ft^){M^)/ ft\M^)?. With our calculation, we can easily 
extrapolate the argument of the denominator, finding the form factor ratio given in Eq. (1.7). 
As a result, Xp = 1.111(94)(34). Reference [20] uses the same sum-rule value A^l^ = 1.24(8) 
when doing the analysis with similar approach, finding the fragmentation-fraction ratio to 
be fs/fd = 0.256(14) (19) (26). We find that 

^ = 0.286(16)stat(21)syst(26)iatt(22)NE, (7.2) 

Id 

where the last two errors (major sources of the theoretical error) are shown explicitly. The 
last error stems from the uncertainty in Me- The result Eq. (7.2) agrees with that of the 
D^TT^ /D^K^ hadronic method, Eq. (7.1), and both agree with LHCfe's determination via 
a method employing semileptonic decays, fs/fd = 0.268(8)stat(t22)syst [16], as well as the 
Particle Data Group's average of LEP and CDF, fjfa = 0.288(24) [15]. 

As a by-product of the calculation, the form-factor ratio in Eq. (1.7) can be combined 
with factorization to estimate the ratio of branching ratios, 

BR{B^ ^ D+K~) ' ^ ^ 

independently of experimental inputs except for quantities like \ Vus\fK/\Vud\fn and lifetimes. 

This work is based on only 4 out of 21 available MILC asqtad ensembles of lattice gauge 
configurations. Further running on ensembles closer to the chiral and continuum limits will 
reduce the length of the extrapolations and, hence, control the growth through extrapolation 
of the statistical error. At the current stage, however, the largest error in Eq. (7.1) remains 
experimental statistics, stemming from the difficulty in reconstructing Df — KKn. 
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Appendix: Staggered chiral perturbation theory for B — >• Div at nonzero recoil 

The material given in this Appendix extends the continuum-QCD xPT for B — )■ Diu [52] 
to staggered fermions. The staggered theory has 16 light pseudoscalar mesons for each meson 
of continuum QCD. Their degeneracy is broken at finite lattice spacing with masses given 
at the leading order by [22] 

Mgg'^ = f^o{mq + mg>) + a^As, (A.l) 



where q, q' are the staggered quarks and /iq is a continuum low-energy constant. a^A= are 
the splittings of the 16 mesons in lattice units, cf. Table II. At this order in an expansion in 
they come in 5 multiplets, labeled P, A, T, V and /, with degeneracies 1, 4, 6, 4 and 1, 
respectively. 

In full (2-1-1) QCD rooted staggered chiral perturbation theory, we have the expression 
for h^-^^: 



,NLO, 



mi 



+ o?5' 



2x' 
V 



+ 



dp* Pit 



Ml 



{Ml-Ml){M, 
Ml - Ml 



f2 



m^-Ml){Ml^-Ml 




If- 

2 



F+ + 



If- 

6 



M2 -M? 

Vv JV 



M2 )(M2 -mm 



(A.2) 
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Similarly for h^^^ {Bs — ?■ Ds), we have 



(s),NLO. 



2x/ 

V 



16 

Ml -Ml 



-Ft + 



Ml -Ml 



m.-Ml){Ml-Ml) ^ {Ml - Ml,^{Ml - Ml^) 



+ 



'Sv '"r?vA^"5v 

Ml - Ml 



-Ft ]+{V->A) 



(A.3) 



where the masses of the fiavor-taste singlet mesons rji and nonsinglet mesons ?7y(^)) Vv(a) 
are given by [51] 

Mrn = I [Ml + 2Ml) , 



Try "I" 5"^ 



4 



V 2 
Z = 



- ( + M^^ + -a^5(, + Z 



(Ml - Mir - la%{Ml - Ml) + A(a^5- 



1/2 



5v Try 

In Eqs. (A.2, A.3), F+ is short for the function F+{w, Mj, A^^V^j)^ defined by 



(A.4) 



F (w, m, x) 



where 



Ii{w,M,x) 



{w + 2)Ii{w, m, x) + {w"^ — l)/2(w, m, x) h{w, m, x) hiui, m, 0) 

(A.5) 

(A.6) 



M 



M'xEi{w) + MV In ( ^ ) ^^(m;) + M'x'Fi{w, x 



and the functions E, G are given by 

E^iw) 



TT 



w + r 

— vr 

(tf + 1)^ 



^3(1/;) = TT, 



with 



r[w] 



2(^2 _ 1) 
1 

2(m;2 _ 1) 
-1, 

1 



[w — r{w)], 

[w"^ + 2 — 3wr(w)], 



ln(rf; + A/ty^ — 1). 



(A.7a) 

(A.7b) 
(A.Tc) 
(A.7d) 

(A.7e) 
(A.7f) 

(A.8) 
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The functions Fi are given by 



where 



F,{w,x) = ^r de Ivrfv^r^-l) -2[/(a)-a]|, (A.9a) 
x'^ Jq l + tysin/6' I \ / J 

FJw,x) = ^ d9- — <^ (Vl-a2-l)H + 

^ ^ x2yo (l + wsin2^)2\ 2 ' 2^1^ 

^"^"V(a)-3al, (A.9b) 



1-a 



(^Vl -a;2 - l) - 2[/(x) - x] I , (A.9c) 



xcos^ 



Vl + U7sin26' 



(A.IO) 



^^""^ i iv/^^^ln[l-2x(x + V^^)] , |x| > li 

and X = A'-^-'/Mj where Mj is the corresponding meson mass. The D(^syD*^^^-^ sphttings are 

^(c) _ ]^4Q g MeV, Ai^"* = 143.9 MeV. When the continuum and chiral hmits are taken, 
the taste sphttings vanish, and the 16 lowest M^s aU tend to the physical pion mass, which 
is around 135 MeV. The extrapolation to the physical pion mass switches from |x| < 1 to 
|x| > 1, requiring both expressions for f{x). 
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